Pearling instability of nanoscale fluid flow confined to a chemical channel 



J. Koplik and T. S. Lo 

Benjamin Levich Institute and Department of Physics, 
City College of the City University of New York, New York, NY 10031 

M. Rauscher and S. Dietrich 
Max- Planck- Institut fur Metallforschung, Heisenbergstr. 3, 70569 Stuttgart, Germany, and 
Institut fur Theoretische und Angewandte Physik, 
Universitdt Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany 
(Dated: February 2, 2008) 

We investigate the flow of a nano-scale incompressible ridge of low- volatility liquid along a "chem- 
ical channel": a long, straight, and completely wetting stripe embedded in a planar substrate, and 
sandwiched between two extended less wetting solid regions. Molecular dynamics simulations, a 
simple long-wavelength approximation, and a full stability analysis based on the Stokes equations 
are used, and give qualitatively consistent results. While thin liquid ridges are stable both statically 
and during flow, a (linear) pearling instability develops if the thickness of the ridge exceeds half 
of the width of the channel. In the flowing case periodic bulges propagate along the channel and 
subsequently merge due to nonlinear effects. However, the ridge does not break up even when the 
flow is unstable, and the qualitative behavior is unchanged even when the fluid can spill over onto 
a partially wetting exterior solid region. 

I. INTRODUCTION 

In recent years substantial efforts have been invested in miniaturizing chemical processes by building microfluidic 
systems. The "lab on a chip" concept integrates a great variety of chemical and physical processes into a single 
device in a similar way as an integrated circuit incorporates many electronic devices into a single chip. These 

microfluidic devices do not only allow for cheap mass production but they can operate with much smaller quantities 
of reactants and reaction products than standard laboratory equipments. Even though most available microfluidic 
devices today have micron sized channels, further miniaturization is leading towards the nano-scale. 

There are two main lines of development for microfluidic systems. The first one encompasses systems with closed 
channels, essentially microfabricated tubes. However, closed channel systems have the disadvantage that they can be 
easily clogged by solute particles such as colloids or large bio-polymers. 

The second type consists of systems which are open with a free liquid-vapor interface, where the fluid is confined 
laterally not by geometrical but by chemical walls. The idea is that the liquid will be guided by lyophilic stripes 
on an otherwise lyophobic substrate |4 til u]- The substrate surfaces can be structured chemically by printing or 
photographic techniques. All these techniques are confined to two dimensions. 

While liquid flow in closed channel systems can be pumped by applying a pressure difference between the inlet and 
the outlet, the liquid- vapor interface of an open system would yield to the pressure (like a very soft rubber tube) and 
droplets would form. Furthermore, upon pumping flow directions in interconnected open channel systems would be 
difficult to predict because small droplets may or may not inflate into larger ones, depending on details of the channel 
geometry and the filling state. However, there are a number of alternative ways to drive a liquid through a chemical 
channel. The most obvious one is to use inertia to drive the flow, e.g., centrifugal forces or gravity, as occurs when 
a liquid film flows down an inclined plane. For films with thicknesses in the nanometer range one needs stronger 
accelerations in order to achieve a decent throughput. These can be realized most conveniently on a rotating disc 
(e.g., a spin coater or a centrifuge). Electric forces are also commonly used in order to drive flows in small channels, 
e.g., by using the electro-osmotic effect which, however, only works for liquids containing ions. Since water is the 
most interesting liquid for biological applications, this is not a serious restriction. Depending on the wall potential, 
a charged layer a few A in thickness forms and an electric field applied parallel to the substrate will set the liquid 
in this layer in motion. While the electro-osmotic effect drives the liquid molecules at the substrate, the Marangoni 
effect drives the molecules at the liquid-vapor interface. The origin of the Marangoni forces are gradients in surface 
tension which arise from temperature gradients or from gradients in surfactant concentrations. The fourth method to 
cause flow in a chemical channel is wicking, i.e., the sucking of a droplet into a chemical channel by capillary forces. 

Here we consider inertial driving of the liquid, which is one of the most promising methods to generate flow in open 
microfluidic devices. Marangoni forces could be induced easily by local heating of the liquid, but for volatile fluids 
like water this will lead to enhanced and undesirable evaporation. Wicking is slow for longer channels and ceases to 
act as a driving force once the channel is filled. In principle, electrical forcing has the advantage that the direction 
of the force can be controlled on very small scales, but the required fields are high and might lead to electrolysis. In 



FIG. 1: A liquid ridge on a chemical wetting stripe of width 2a embedded in a solid substrate, which is non-wetting. Gravity or 
inertial forces act on the liquid in the direction parallel to the z-axis. The liquid-vapor interface is parameterized in cylindrical 
coordinates as r = R(<f),z,t), and the local height h(z,t) — R(n/2, z,t). 



contrast, inertial driving can be realized easily by centrifugal forces and does not affect the liquid sample. 

Our focus in this paper is the stability and robustness of one-dimensional channel flow (i.e., flow in a straight 
channel) with inertial forcing aligned with the channel direction. The flow geometry is sketched in Fig. ^ In the 
absence of forcing, a chemical channel filled homogeneously by a non-volatile liquid of fixed volume is unstable with 
respect to droplet formation when the contact angle with the triple-line pinned at the channel edge exceeds 90° 0, • 
This pearling instability is surface tension driven and similar to the Rayleigh-Plateau instability. The resulting drop 
shapes have been studied extensively in Ref. . Of particular interest is the question of spilling onto the lyophobic 
substrate, which can lead to cross-talk between neighboring channels. In the force-free case this happens for large 
droplets on chemical stripes on a partially wetting embedding substrate. 

In Sec. [n] we describe the formulation of the molecular dynamics (MD) simulations for flow on a chemical channel 
and discuss the special case of unidirectional flow quantitatively. In Sec. IIIII we discuss the appearance of pearling 
instabilities in channel flows and, in particular, we study the influence of the liquid ridge height (effectively, the 
channel depth) and of the wetting properties of the surrounding substrate on the stability of the liquid ridge. The 
full linear stability analysis of a Stokes flow in Sec. IIVI as well as the hydrodynamic long- wavelength approximation 
developed in Sec. are in qualitative agreement with the MD simulations. We discuss the results in Sec. I VII 



II. MOLECULAR DYNAMICS SIMULATIONS 



The MD simulations carried out here employ standard techniques [Tol 111) and a simple molecular model for a 
non-polar liquid in contact with a thermally-fluctuating solid [I2>U3> based on atoms interacting via an adjustable 
Lennard- Jones potential 
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The interaction is cut off at r c — 2.5 a, and shifted by a linear term so that the force vanishes smoothly there. In this 
section we employ so-called "MD units" derived from this potential, with the atomic core size a as the unit of length, 
the potential well depth e as the unit of energy, and r = o^m/e) 1 / 2 as the unit of time t, where m is the mass of 
the liquid atoms. We simulate a generic non-polar liquid rather than any particular laboratory material, but typical 
numerical values are a w 0.3 nm, e/fcs ~ 140K, and r ps 2 ps. A monatomic Lennard-Jones fluid exhibits a high 
vapor pressure and a very broad liquid-vapor interface. These are inconvenient properties in the present context and 
are ameliorated by using a FENE potential 



V F (r) = -~k F r 2 ln( I - - 
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FIG. 2: Snapshot of a typical initial molecular configuration. The left and right views are along and top down on the stripe, 
respectively. The liquid molecules are depicted as the three "bonds" joining the four atoms, while substrate atoms with wetting 
(cls = 1) and non-wetting (cls = 0) properties are light green and dark red dots, respectively. The average spacing between 
solid atoms is 0.78a. The liquid ridge has a circular cross section (as indicated by the dashed circle) and a height ho (at which 
the density attains half its value in the center), which depends on the liquid volume and the stripe width 2a. 

to group the atoms into freely-jointed chain molecules which are four atoms in length. The numerical values of the 
constants in Vf are kp = 30e/cr 2 and ro = 1.5c, following Ref. [14j . This leads to a significant reduction of the vapor 
pressure and thus yields a liquid of low volatility. The solid atoms are tethered to two layers of fee lattice sites {roi} 
with lattice constant 1.55 a using a stiff spring 

1 2 

Vs(ri) = -k s Oi -r 0i ) i = l,...,n s (3) 

with spring constant ks = lOOe/er 2 , allowing the liquid and solid to exchange energy and momentum. The substrate 
surface has a (100) crystallographic orientation. The solid atom mass is chosen as 100 m, so as to have approximately 
the same oscillation frequency as the LJ interactions and a common time step in the numerical integration. During 
the simulations, a Nose-Hoover thermostat maintains a constant temperature l.Oe/ks- Under these conditions, the 
fluid condenses into a liquid of number density 0.79 cr -3 , in equilibrium with a modest amount of vapor. 

The coefficient in Eq. Q is used to vary the strength of the attractive interaction between atomic species i and 
j. The intra- liquid coefficient has the (standard) value cll — 1-0. With the other modeling choices and parameter 
values given above, the liquid wets the solid completely, partially, and not at all for the choices cls = 1.0, 0.75, and 
0.0, respectively. The corresponding contact angles for a drop on a homogeneous solid with these values of cls are 
0°, approximately 90°, and 180°, respectively. In the latter case and in the absence of gravity, a drop actually drifts 
off the surface. The surface tension of the liquid is found from a standard simulation of a slab of liquid coexisting 
with vapor to be 7 = 0.46 e/<7 2 , and the liquid viscosity is obtained from a separate simulation of Couette flow as 
i] = 3.60 to/ (err). 

A snapshot of the basic simulated configuration is shown in Fig. [21 i.e., the individual atoms in a liquid ridge on a 
wetting stripe for a case where the exterior of the channel is completely non-wetting. The left figure shows an end-on 
view (parallel to the stripe) and the right figure a section of the top view. Note that the interaction does confine the 
liquid to the wetting region, and that the interface strongly fluctuates. In this and all other simulations, we apply 
periodic boundary conditions in the plane of the substrate in order to avoid end effects. To construct the atomic 
configurations, we place atoms on the sites of a cylindrical section of a fee lattice and allow them to move under 
the action of the L J and FENE potentials given above until the atomic positions become disordered and the internal 
energy stabilizes. This configuration does not correspond to global thermodynamic equilibrium because, as we shall 
see, a uniformly shaped liquid ridge is hydrodynamically unstable; however, this configuration can be described in 
terms of a local equilibrium. The resulting cylinder of liquid is translated to lie just above and parallel to the wetting 
stripe of a solid substrate which has been independently allowed to reach local equilibrium, so that the liquid is 
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attracted towards the substrate and adjusts itself to form an irregular circular ridge with height ho as shown in Fig. [21 
During the course of the equilibration process, lasting typically several hundred r, the liquid surface fluctuates at a 
5% level in an apparently random manner. Systematic changes in the shape of the ridge occur only on longer time 
scales. 

In order to simulate liquid flow, we apply a uniform force mge z to each liquid atom, where the axis of the 
stripe is the z-direction (see Fig. The results presented in this paper involve values g = 0.01-0.05 cr/r 2 , and the 
corresponding Reynolds numbers, based on the observed mean flow velocities and the height of the liquid ridge, are 
Rc = 0.1-10.0. Smaller values of g are perhaps more realistic but require longer computations for the flow and for any 
changes in interfacial shape to develop. In any case, these low-Reynolds number simulations are in the linear regime 
where phenomena simply scale in time with l/g. At higher g, the liquid is often driven off the substrate. 

In order to relate the MD simulation results to macroscopic hydrodynamics, we consider the special case in which 
the liquid occupies just half of a cylinder with a 90° contact angle, with initial height h$ — a. As discussed below in 
Sects. IIVI and IVl this case is neutrally stable with respect to changes in interface shape and the shape of the ridge 
does not change significantly over the time of the simulation. Furthermore, this case may be directly compared (over 
moderate time intervals) with solutions of the Navier-Stokes (NS) equations for unidirectional flow parallel to the 
stripe axis. Furthermore, the unidirectional NS flow field may be computed in an almost closed form by a Fourier 
series method (see Sec. IIVI for details). In Fig. [3] we show the liquid density, the axial velocity field, and the r-z 
component of the shear stress, as obtained from the analytic solution and from the MD simulations, at the same 
resolution. (The explicit analytic form of the shear stress can be obtained from the formulas in Sec. IIVI and the MD 
counterpart is the Irving-Kirkwood expression ^3-) ^ n order to obtain these fields from the simulations, the flow 
domain in the x-y plane was divided into square bins of size a x er, and the appropriate quantity in each bin was 
averaged both over time (2000 t) and the length of the system in z-direction (517 a). A larger averaging time interval 
was impractical due to the eventual development of systematic interfacial fluctuations, while smaller sampling bins 
lead to excessive noise in the stress (and wall- induced oscillations in the liquid density near the solid [12|). The 
jaggedness in the NS figures results from the relatively small number of bins (31 x 16) occupied by the liquid, while 
that in the MD figures incorporates thermal fluctuations as well. 

The first point to note is that the MD liquid-vapor interface is rather diffuse, corresponding to a width of about 7 
sampling bins. This feature is not easily accounted for in the NS calculation, so that close agreement should not be 
expected. Secondly, we see that while the shear stress values are in qualitative agreement both in shape and magnitude 
and the velocity fields have the same shape, the velocity values differ by almost a factor of 2. This comparison depends 
critically on the value chosen for the liquid ridge height in the NS calculations: a larger height increases the velocity 
values and decreases the stress, so that a different choice could make both agree to 50%. Here, in computing the NS 
flow field, we selected the numerical value of the height ho based on the reasonable but arbitrary definition as the 
point where the MD density dropped to half the value at the center of the ridge (see Figs. ^ and [2]). The fact that 
the MD velocities are larger may result from the fact the free surface region of the ridge exhibits a transition between 
liquid and gaseous phases and time-averaged velocities in the latter are larger. In any case, the degree of agreement 
between MD and NS hydrodynamics is worse than in similar comparisons where the liquid is completely confined 
by walls 01 > which we attribute to the diffuse top boundary of the liquid. A final remark is that the (thermal) 
fluctuations in the MD results are largest for the stress, moderate for the velocity, and almost invisible in the density. 
The reason is that the stress computation involves the intermolecular force, which is a rapidly varying function of 
atomic separation, while the velocities result from an integration over the force, which provides some smoothing. The 
density follows from the atomic positions, which are a further integral of the velocity, and is thus still smoother. 

III. PEARLING INSTABILITIES 

Turning now to the systematic study of pearling instabilities, we have conducted a number of simulations in which we 
explore the effects of several operating parameters. For a fixed number of liquid atoms (128,000, or 32,000 molecules) 
we vary 

• filling: the number of liquid particles per stripe length on the wetting stripe. In practice, we choose different 
widths for the solid wetting stripe, and since the liquid is effectively pinned to the edges of the stripe, it adjusts 
its height and contact angle such that a circular cross section emerges. 

• wetting: the region exterior to the completely wetting stripe can be either completely non- wetting (i.e., drying) 
or partially wetting. In the former case the liquid is confined to the stripe and does not contact the exterior, 
but in the latter case it can spill over the edges. 

• forcing: the value of the acceleration g. 
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FIG. 3: Comparison of the density, axial velocity and r-z shear stress fields in a liquid ridge with initial height ho = a = 16a 
for unidirectional flow with g = 0.01 (see Fig. Qfor the coordinate definitions). In each frame of the figure the base is the x-y 
plane, with the substrate (not shown) to the left. The height of the surface represents the value of the field depicted there, 
and the grid spacing is a. In the density plots, the plateau height, representing the central liquid density, is 0.79 a~ 3 . In the 
Navier-Stokes (NS) velocity and stress plots, the maximum values are 0.225 ct/t and 0.107 m/(o-r 2 ), respectively, and for each 
field the same scales and units are used for NS and MD. The MD fields are set to zero in bins (of size a x er) where the density 
is below 0.01 a~ 3 . The exterior of the substrate is nonwetting (cls = 0) whereas the stripe is wetting (cls = 1). 
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FIG. 4: Head on view along the 2-axis of the initial configuration of the cases simulated. From left to right: non-wetting 
(cls = 0) exterior and narrow wetting (cls = 1) stripe with ho/a = 4.0; non-wetting exterior and wide wetting stripe with 
ho/a = 2.2; partially wetting (cls = 0.75) exterior and narrow wetting stripe with ho/a = 3.1; partially wetting exterior and 
wide wetting stripe with ho/a = 1.7. The number of molecules per channel length is the same in all cases. 



• length: simulations have been carried out using one-half and one-quarter of the stripe length and, correspond- 
ingly, of the number of atoms, but with identical cross-sectional shapes. This allows one to examine the 
wavelength selection if pearling develops. 

At early times before any instability has appeared, Fig. 0] gives head-on snapshots of the cross sectional shape of 
the liquid for the two choices of filling and two choices of wetting considered in this paper. In the two cases of a 
completely non- wetting exterior solid shown in Fig0J the key geometric parameters for the Stokes equation analysis 
below are the wetting stripe width 2 a and the initial liquid height ho- The case of a narrow stripe corresponds to a 
width 2 a = 10.6a and a height-to-half-width ratio ho/a — 4.0, while the case of a wide stripe corresponds to a width 
2a = 17.1a with ho/a — 2.2. The maximum system length considered here is 547 cr with 128,000 atoms. 

First we consider the case of no forcing. If a static liquid ridge of uniform cross-sectional shape with contact lines 
pinned at the chemical steps (ridges with a mobile contact line are always unstable Q)is sufficiently long and thick, 
it is expected to be unstable to the formation of bulges, as shown by various authors [1,01 EE]- This instability is 
driven by surface tension, and is closely related to the well-known Rayleigh-Plateau instability of a liquid cylinder, 
which reduces its surface area under long-wavelength perturbations of its radius, eventually breaking up into spheres. 
Here, the qualification "thick" means that the contact angle must be larger than 90° or equivalently, because the 
liquid cross section in equilibrium must be an arc of a circle, the initial height must be greater than half the width of 
the wetting region, i.e., ho > a. 

This static instability appears naturally in MD simulations. In Fig. [3] we show an example where the initial height 
was 2.0 times the stripe width, and the system had 64,000 atoms and half the maximum length studied by us. In 
this and in subsequent figures, instead of individual atoms we plot the median liquid interface, defined as the surface 
on which the liquid density averaged over a 50r interval falls to half of the value at the center. A second simulation 
with initial height ho = a was stable over the same time interval, showing continued random fluctuations in shape 
but no systematic evolution. Because the liquid is strongly attracted to the wetting region of the substrate, the 
pearls (with a characteristic wavelength) arising from the instability remain connected by thin liquid films over the 
duration of the simulation. Neither the simulations, nor the linear stability analyses below, can reliably predict the 
ultimate asymptotic state and we cannot determine whether the liquid pearls remain connected or not in the true 
asymptotic sense t — > oo, i.e., in thermal equilibrium. Note that the Surface Evolver software used in |5ll^.ll5j| does not 
incorporate the substrate potential as such and cannot address the issue of adsorbed films of molecular thicknesses. 
This could be resolved by density functional techniques combined with a constraint for the liquid volume; we regard 
it as likely that the pearls remain connected by thin liquid bridges. In the simpler case of the pure Rayleigh-Plateau 
instability, i.e., in the absence of a substrate, MD simulations similar to those described here are in rough agreement 
with theoretical predictions of the breakup time; see Refs. \lh\ \lj\ . 

The behavior of the same system when driven with an acceleration g = 0.01 is shown in Fig. for the same 
initial shape. Also in this case pearls form, with the same initial wavelength as in the unforced case, and then 
propagate downstream with the moving liquid. While initially the deviations from the uniform thickness of the ridge 
are periodic and have the same velocity, as they develop to finite amplitude they presumably interact nonlinearly 
and acquire different velocities, and eventually merge into pairs until one moving pearl survives. If the length of the 
system is half as long, initially only two pearls are present, which eventually merge into one; with one quarter of the 
length, one pearl appears and it simply grows and propagates. If the simulation is repeated with g increased to 0.025, 
the same instability scenario is present (see Fig. [7J| but the structure develops more rapidly, the pearls propagate at 
a higher velocity, and mergers occur sooner. At still higher values of acceleration, in the inertial regime with Rc > 1, 
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FIG. 5: Instability of a static liquid ridge with ho /a = 4.0 in the nonwetting- narrow case (see Fig.^J. The sequence runs from 
left to right, showing the median interface position at the times indicated. The stripe length is 273.6a and the vertical axis at 
the rear of the figure has a height of 51.3a. 




FIG. 6: Time evolution of the instability of a driven liquid ridge with the same aspect ratio and substrate properties as in 
Fig. |5] (i.e., ho/a — 4.0) but with twice the system size driven at g = 0.01, showing propagating pearls and their merger. The 
forcing is along the channel from back to front, and the vertical axis at the rear of the figure has a height of 51.3<r. 

the liquid tends to assume very irregular shapes and may even lose contact with the wetting stripe. 

Somewhat surprisingly, the case that the embedding substrate exhibits partial wetting shows an identical qualitative 
behavior in terms of pearl formation, motion and merger, with only numerical differences in growth rate, wavelength, 
and velocity. Generally speaking, in these cases the motion is slower because the liquid spills over laterally and tends 
to be closer to the solid substrate, because the exterior part of the substrate is more attractive than in the previous 
case. The same number of pearls appear initially but require a longer time to develop and to subsequently merge. 
An example is shown in Fig. [H] for the case of a wide stripe of width 2 a = 17.1a at g = 0.01, and the principal 
difference from the corresponding results for a non-wetting exterior is the longer time scale. The corresponding top 




FIG. 7: The same system as in Fig. [5] but driven at a higher acceleration g = 0.025. Note the difference in time scale for the 
merging of the pearls. The vertical axis at the rear of the figure has a height of 51.3c 
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TABLE I: Comparison of MD simulation, long wavelength (LW) approximation, and Navier-Stokes (NS) equations stability 
analysis, for four cases with a non- wetting exterior substrate. Two stripe widths (wide, ho /a = 2.2, and narrow, ho/ a = 4) and 
two accelerations g were considered. T is a rough estimate of the time at which pearls start to become observable in the MD 
simulations, which is compared to Wr 1 , the inverse of the real part of the instability growth rate. A* and v are the distance 
between the centers of neighboring pearls and their velocity, respectively. 




FIG. 8: Spatio-temporal effects due to a partially- wetting (cls = 0.75) exterior solid for a wide wetting stripe (ho /a = 2.1) at 
g = 0.01. Left: perspective view, right: top view. The vertical axis at the rear of the left figure has a height of 51.3a, and the 
contours in the right figure are lines of constant height. The scale bars indicate the width 2 a of the wetting stripe. 

views demonstrate that the lateral contact line formed by the liquid on the exterior solid tends to move outward near 
a pearl, but it never deviates very far from the edge of the stripe. 

A summary of the numerical results for the onset time of the instability, its initial wavelength, and the pearl velocity 
for the cases with a non- wetting exterior substrate is given in Table din the rows labeled MD, while the other entries 
are obtained from the stability analyses given below. In the table, the "MD" numbers assigned to the instability 
are rough estimates obtained by visual inspection. The onset time T is the time elapsed after the application of the 
forcing when a recognizable periodic disturbance (amplitude roughly 10% of ho) appears in a side view of the liquid 
ridge. The wavelength is the average spacing between crests in the pattern, and the velocity is obtained from the 
displacement of the crests of the pearls at successive times. The numerical entries for the two analytic calculations 
are the values for the most unstable mode in each case. The situations involving a partially wetting exterior substrate 
present a rather more difficult stability problem, because in that case the lateral sides of the liquid ridge are mobile 
instead of being pinned at the stripe edges. We have not undertaken this stability analysis. 

Of course, it would be preferable to compare MD and continuum results in terms of the time dependence of individual 
Fourier components of the deviation from a uniform axial shape. But in practice two difficulties prevent this analysis. 
First, on the atomic scale a, the simulated liquid ridge exhibits a diffuse interfacial region, so that a prescription is 
needed to define a unique liquid-vapor interface. Second, independent of how the interface is denned, it undergoes 
capillary wave-like thermal fluctuations which can obscure the relevant signal. The procedure we tried was to choose 
equally spaced time intervals, to divide the liquid ridge into slices of constant-thickness along the flow axis, and at 
each time to count the number of atoms in each bin. Assuming that the number of atoms in each bin is proportional 
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to the cross-sectional area at that position z, and that locally the interface is an arc of a circle, geometry gives the 
local height h(z,t) (see, c.f., Eq. 1|26[) '). The power spectrum of h is obtained by a fast Fourier transform, from which 
the time evolution of individual Fourier modes follows. If these modes had displayed an exponential growth at early 
times, a precise comparison with a stability analysis would be possible. However, we did not observe any systematic 
behavior of the modal growth rates: some increase with time but not exponentially over any appreciable time interval, 
while others simply vary in an irregular manner due to thermal fluctuations. In contrast, a recent MD calculation of 
the Rayleigh- Taylor instability in a two-liquid system was able to relate the interfacial shape to the growth of distinct 
Fourier modes 18]. The presence of the second liquid sharpens the interface and suppresses thermal fluctuations, 
allowing one to make closer contact with the spectral analysis. 

Finally, we consider the effects of the pearling instability on liquid transport. One effect of pearl formation is to 
cause a net displacement of fluid molecules away from the substrate where they are partially bound by the substrate 
atoms providing a wetting condition, so that the mean velocity and the flux should increase as the instability grows. 
From a continuum point of view, the equivalent statement is that pearling moves the liquid further away from that 
part of the boundary on which the no-slip condition applies. In fact, as a function of time we observe a roughly linear 
increase of the mean axial velocity (w) (averaged over the whole liquid ridge) as shown for one example in Fig. 
The velocities with which the pearls move are comparable to the mean flow velocities in each case. In view of the 
approximate way in which the pearl velocities are determined, i.e., by monitoring the position of the crests at regular 
intervals in time, it is difficult to obtain a more precise statement. At longer times, the mean velocity oscillates around 
the maximum value in Fig. reflecting shape adjustments of the surviving pearl. Since in this case Re = 0.1, the 
roughly linear increase of (w) with time can be attributed to shape changes rather than to inertia. 

2.5 
2 

J 1 
0.5 


500 1000 1500 2000 

t/X 

FIG. 9: Mean axial velocity (w) vs. time for the case of a narrow wetting stripe embedded into a non-wetting substrate 
(ho/a — 4) as shown in Fig. [(Jdriven at g — 0.01. 




IV. STABILITY ANALYSIS 



In this section, we analyze the linear stability of the uniform flow state on a wetting stripe in the framework of 
continuum fluid mechanics. For simplicity we consider stripes of infinite length in order to avoid the problem of 
wavelength selection. Assuming the liquid is Newtonian and incompressible, its flow velocity u = ue r + ve^ + we z for 
the geometry shown in Fig. ^ is governed by the Navier-Stokes equation 

— + u-Vu = — Vp+^V 2 u + .ge, 4 
ot p p 

and the incompressibility condition V • u = 0, where p and p are the mass density and the pressure, respectively. In 
addition, the deformation of the free surface R((f),z,t) is described by the kinematic condition 

OR v dR OR 

-dt + Rd6 +W ^ =U > (5) 



which expresses the fact that (for a non-volatile liquid) a fluid particle at the interface moves inwards or outwards 
with the interfacial velocity. At the impermeable solid, u satisfies the no-slip boundary condition (which is consistent 
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with the MD simulations). At the free surface, assuming that the liquid is non- volatile, the appropriate boundary 
conditions are the balance of the normal component of the stress n-(T-n) = — 7/t, and the condition of zero tangential 
force, t-(T-n) = 0. Here k is the local mean curvature of the liquid-vapor interface, T = — pi + ?7[(Vu) + (Vu) T ] is the 
stress tensor, n is the unit normal vector pointing towards the vapor phase, and t with n • t = is any unit tangent 
vector on the surface. Since we are considering a non-volatile liquid, we can assume that the density, the viscosity, 
and the pressure of the vapor phase are negligible. 

Referring to Fig. ^ we assume that surface tension and substrate potential act - for a sufficiently large liquid 
volume - to pin the liquid at the edges of the wetting stripe, and to fix the initial liquid-vapor interface to be a 
(constant-curvature) section of a circular cylinder. The contact angle at the edge of the liquid can then take on any 
value between that corresponding to the wetting region on the stripe (0°) and that to the non- wetting one on the 
exterior (180°). For any given constant initial height h and acceleration g, Eq. (0} has a steady-state unidirectional 
flow solution u = wq(t, <fi) e z which satisfies the two-dimensional Poisson equation 



1 d_ fdwo\ 1 d 2 w _ _pg 



r 



dr \ dr J r 2 d<j) 2 r\ 



(6) 



The boundary conditions for u reduce to wq — at the substrate (|x| < a, y = 0) and dwo/dn — at the free surface 
(r = R, y > 0). The axially constrained shape of the initial liquid- vapor interface can be written as 



R(<f>, z, t) -> Ro(<j>) = a {b„ sin <f> + y 1 + b 2 sin 2 <f>) (7) 

with 

-Kt-£)- 

The pressure inside the liquid is po = j/R c where R c = (h 2 + a 2 )/ (2 ho) is the radius of the truncated cylinder. For 
general ho/ a, we obtain the base state Wq by solving Eq. 10 numerically, using a finite difference discretization. In the 
special case in which the liquid occupies exactly half of a circular cylinder (i.e., h /a = 1), one can obtain an explicit 
Fourier series solution with the separation ansatz wo (r, </)) = — pgr 2 /(Ar]) + Bo + Y^=i r ™ [A n sm(n(f))+B n cos(n (/>)]. 
The boundary conditions and the symmetry of w (r, </>) determine the coefficients A n and B n , and we obtain 



wo{r, <p) = 

V 



1 \ - S (T 



L / r \ 2 

4 U (C ° S2 ^ 



^— ' tt n 2 (n 2 — 4) \a 

n=l,3,5... V ' 



— sinnc 



(9) 



which has been discussed in Sec. HP and is plotted as the NS velocity in Fig. 01 

If the liquid starts from rest, the time scale for the steady state to develop is controlled by the diffusion of vorticity, 
which happens in the cases simulated above on the time scale td — pa 2 /r) w 6 — 16 r, which is much smaller than 
the characteristic time for the long wavelength instability to appear. It is reasonable, therefore, to study the stability 
of the flowing liquid ridge by investigating a small perturbation about the basic unidirectional steady state. Hence, 
we write u = u'e r + v'e^ + (w + w')e z , p = p + p' and R = Ro((j>) + R'{4>, z, t) and linearize Eqs. Q, 0, and all 
the necessary boundary conditions about the basic state. This leads to a set of linear partial differential equations 
(PDE) for the small quantities u', v', w' , p' and R' . The results are fl9j : 

m = A d A + a (Vv - - - -—) - wo— (10) 

dt p dr p \ r 2 r 2 dcf> J dz 



dv' 1 dp' i r] ( 2 1 v' t 2 du'\ dv' 



for the velocity and 



dt pr d(j) p \ r 2 r 2 dip) ° dz ^ ^ 

dw' 1 dp' 77 „, , dw' , „, 

— = —JL + lTpw'-wz— 12 
dt p dz p dz 



dR' , v' dR a 8Bf 

= u _ ~5~~^T ~ w o^r~ (13) 



dt Ro d4> dz 

for the kinematic condition in Eq. |JH}. At any time t, p' satisfies the pressure equation 



vy = -2 P (^ + i%^) , 

1 dr dz r dip dz 
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FIG. 10: Real (left) and imaginary (right) parts of the growth rate u)(k, 1) (in units of td = pa 2 /rj) of the dominant growing 
instability modes as function of the wavenumber k of the perturbation and for n = 1 (see Eq. 1221 1 for a liquid ridge with 
ho/a — 2.0. The wavelength of the fastest growing mode and the corresponding growth rate increase with g. The dashed line 
marks the loci of these maxima as g varies from g = to 0.1 from bottom to top. The accelerations g are measured in units of 
<jt~ 2 . Note that cjR(k > fc max (g)) < 0, which implies that perturbations with shorter wavelengths are linearly stable. Stronger 
forcing leads to faster growth of the most unstable perturbation but in the same time leads to a stabilization of perturbations 
with shorter wavelengths. 



which is obtained by taking the divergence of Eq. and using the incompressibility condition. 

The velocity boundary conditions, upon evaluation at the liquid- vapor interface position r = i?o(^>) + R'( 
and after linearization, are: 



du' 
dz 



dw' 
dr 



Ro4> ( <V 
Ro \dz 



from V(T-n) = 0; 



1 dwo 



1 dw n 

~R~ lty ) Rl d<p 



Ra 



2 -^W R i 
Ro 



d 2 w Q _ Rp^ d 2 w 
dr 2 i?,g drd(f> 



R' = 



(15) 



1 - 



R o<t> 
Rl 



ch/_ 
dr 



1 

Ro 



du 1 



2R, 



00 



Ro 



} du' 
dr 



dw' 
dz 



1 

Ro 



1 



R l<t> 
Rq 



dw 



R' = 



(16) 



from t^-(T-n) = 0; and 



du 1 
dr 



Ro 



dv' 



dw 1 
1h 



(17) 



from V u = 0. The no-slip boundary condition at the substrate requires u' = v' = w' = at <fi = and 
Applying the normal stress condition at the liquid-vapor interface r — Ro((p) + R'{4>, z, t) gives 



<f> = IT. 



Rj,. 



(l + «oV^o 2 ) 1/2 



Rl(l + Rl JEZ)W 



Ro ( 



Ro4>4> 
Rn 



4^ 
Rl 



1 + Rl^/Rl 

2r, du' R 04> (dv' 

dr R \ dr 



RL 



1 ~ 2R l<t>/ ^0 



1 du' 



Ro Ro 



R 0<t> 
R 3 



dv' 



2R o<t> 
R 2 



R 



000 



Ro 



R' 



after linearization. In Eqs. I|15|) . (|16|l . and l]18p. the subscript 4> denotes differentiation with respect to 
r = R (4>). 



(18) 
and 
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We have solved the linearized equations using the pressure Poisson equation (PPE) formulation |20ll2l|. which re- 
quires boundary conditions for the pressure p' at the substrate. These conditions can be obtained from the projections 
of Eq. (@J along n and t, and from the incompressibility condition. To linear order, they become 

I dp' / 2 , 2 du'\ 

dp' _ / a , 2 ft,' 



and 

1 d(ru') 1 ft,' 3u>' 
r ft r ft/> ft 



(21) 



at = or 7r. 

We represent the small quantities in terms of Fourier series expansions, and finally derive a set of evolution equations 
for the Fourier components. In this way, we reduce the three-dimensional problem to a set of two-dimensional linear 
PDEs for each Fourier mode, characterized by a wavenumber k. We consider a small initial sinusoidal perturbation of 
the liquid-vapor interface of the form R'(<j>,z,t = 0) = e cos kz sin <f> for |e| <C Ro, and follow its evolution by solving 
the initial value problem numerically. The solution is expanded in a Fourier series of the form 

OO r. 

R'{(j>,z,t) = ^2 dk A nk(t)e- ikz sinn0 , (22) 

n=l 

and the real and imaginary parts of the dominant growth rate Lo[k, n) = LUu(k,n) + iuji(k,n) are extracted from 
the time-dependent amplitude A n k(t) oc e u ^ k ' n ' t . Modes with n > 1 are stable. The functions uj(k, 1) for the case 
ho /a = 2 at different values of g are shown in Fig. ^3 It is evident from Fig. ^] that both the amplitude growth 
rate lur and the phase velocity uij/k increase as g increases, which is consistent with the MD simulations. The slight 
change in wavelength of the maximally unstable mode between g = 0.01 and g — 0.025, which are the values used in 
the MD simulations, could not be detected in the simulations for the reasons discussed in Sec. [H] We have carried 
out the same calculation also for the case ho /a = 4 and the results are qualitatively the same. For the case ho/a = 1, 
all Fourier modes are stable. 



V. LONG- WAVELENGTH APPROXIMATION 



Since the initial liquid configuration is a long ridge with a uniform cross section, and the structures which develop 
have, at least initially, characteristic length scales much larger than the width of the ridge, it is natural to consider 
a long-wavelength approximation to the Navier-Stokes equation. Referring to Fig. ^ incompressibility implies the 
continuity equation 



OA _ _dQ 
dt dz 



(23) 



where A is the local cross-sectional area 

i r 



A(z,t) = ^ I d4>R 2 (^z,t) (24) 



and Q is the volumetric flow rate given by 



Q(z,t)= I d<t> I drrw(r,(/),z,t) (25) 



o Jo 



with w(r, <p, z,t) being the axial component of the flow velocity. 

If we assume that the local cross section always has a circlular boundary, the free surface can be described by 
Eq. Q with b — b(z, t) expressed in terms of the local height h — h(z, t) according to the same functional relation as 
between the initial value bo and the initial height ho, Eq. JSJ. In this case, the integral in Eq. I|24|) yields 



A = a A 



(1 + b 2 ) +arctanb) , (26) 
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h/a h/a 

FIG. 11: The mobility fi(h/a) as defined in Eq. (13011 in units of (h/a) 4 (left) and its derivative in units of (h/a) 3 (right) from 
a numerical solution of the Stokes equation. The dashed lines correspond to the thin film limit h/a — > 0. The value of /i(l) 
is known exactly. These units are chosen as to make the difference between the numerical solutions and the thin film limit 
particularly visible. Moreover, it is the combination a 4 fi(h/a) which enters into Eq. l-'iDI for the flux Q. 



so that 



OA / a 2 



I I) ( — + arctan b 



Oh 
~0t 



(27) 



Instead of performing a systematic long wavelength expansion we further assume that w(r, (p,z,t) can be approxi- 
mated by the steady state axial flow field of a liquid ridge with uniform cross-sectional shape given by the local profile 
R(4>, z, t). The flow is driven by the effective forcing (g — which means that w(r, (j>, z, t) is given by the solution 



of Eq. © with m replaced by ( 



idp\ 

r\ dz>' 



This allows for the ansatz 



w(r, cf>, z, t) = wi (r, 4>; z,t) 



pg 
n 



1 dp 

r\ dz 



(28) 



due to the linearity of the Poisson equation. Here w% has the dimension of an area and is proportional to the base 
state velocity wq of the previous section because it obeys the same boundary conditions. Note that although w\ is 
proportional to the velocity in a uniform infinite ridge, a function of r and 0, it depends implicitly on z and t because 
of the boundary conditions imposed at r = R(4>, z, t). The Laplace pressure p in Eq. (|28|l is approximately given by 



2h 



h 2 



d 2 h 

dz 2 



(29) 



where the first term 2h/ (ft 2 + a 2 ) is the inverse radius of curvature of the local circular cross-section set by the local 
height h, and the second is the (leading order approximation to) the curvature in the orthogonal plane. Within this 
approximation the flux is 



Q 



rRoW 

d(f> / drrwi(r,(f>;z,t) 
ii Jo 



Pg 



1 dp 
r\ dz 



a ,h 
— fli- 
rt 



a 



pg 



dp 
dz 



(30) 



Since wi is proportional to Wq, the solution for w\ can be inferred from Sec. IIVI and in general it is obtained 
numerically. The dimcnsionlcss "mobility" fi(h/a) (i.e., the proportionality factor between the flux-driving pressure 
gradient and the flow rate), defined via Eq. I|3U|) . is obtained by numerical integration of the solution of Eq. In 
Sec.|n]we considered the special case h/a = 1, for which Eq. provides the analytical solution, which in turns gives 
/i(l) = | {— [6 + 7£(3)] — 7r} = 0.362. The thin film limit h/a — >• is an interesting special case as well. In this limit, 
we have a Poiseuille velocity profile locally and therefore Wi = Wi(x,y) = £(x)y — \y 2 where l(x) = h [l — (x/a) 2 ~\ is 



the local film height as a function of the lateral position x. This leads to fi(h/a — ► 0) = (h/a) . 
the mobility /i and its derivative as functions of h/a together with their limiting behaviors for h/a 
With Eqs. H27fl and (|30|l . the continuity equation (|23|l can be re- written as 



1 + 



arctan b 



dh 
If 



,h. d 2 p 1 ,,h. 
a' dz z a a' 



Op 

dz 



P9--5-)-W 



dh 
dz 



Figure ^2 shows 
0. 



(31) 
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Denoting the local height by h = ho 

„2 



Sh(z, t) and linearizing Eq. I|31|) around ho gives 



1 



h 2 
"o 



l + b 



arctan bo 



dSh 



7 a 



r\ v a ' 



2 (h 2 - a 2 ) d 2 5h d 4 5h 



(h 2 + a 2 ) 2 dz 2 dz* 



(32) 



pga 2 ,,h ,d5h 



with h 



o 



(ho/a 



rj a ' dz 

a/ho)/2. Specializing now to a sinusoidal perturbation 5h ~ 

From this, one can identify the real and imaginary parts of the growth rate u>: 
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1) 
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(33) 

(34) 
(35) 



Graphs of ujr for the two cases simulated by MD are shown in Fig. ^| Equation (|34|) shows that a long wavelength 
perturbation to a liquid ridge of uniform height is unstable for ho > a (i.e., for contact angles larger than 90°) 
and is stable otherwise. Figure 1121 indicates that for thick ridges the instability grows more rapidly but only for 
longer wavelengths. Equation l|35|l shows that the phase velocity loijk is proportional to the acceleration g, which is 
consistent with the observation in MD that stronger forcing leads to faster moving pearls. However, within this simple 
approach, the amplitude growth rate lor is independent of g opposite to what is observed in the MD simulations and 
predicted by the stability analysis in Sec. II VI Numerical values for the growth rate, the critical wavelength, and the 
phase velocity are compared with the simulation results in Table I, showing at least qualitative agreement. 
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FIG. 12: The real part of the growth rate of unstable modes for two values of ho /a in the long- wavelength approximation (see 
Eq. 13410 . In contrast to the full linear stability analysis in Sec. II VI the real part does not depend on g. lor is negative for 
ho j a < 1. 



VI. SUMMARY AND DISCUSSION 



Motivated by possible applications in the field of open microfluidic systems we have studied the flow of a non- 
volatile liquid ridge along a stripe- like chemical channel, to which the liquid is confined by wettability (see Fig. ^| . In 
this context the most salient question concerns stability, i.e., the degree to which the liquid remains on the channel 
and the effects of any changes in the shape of the liquid-vapor interface. Our basic tool to calculate the flow has 
been large scale molecular dynamics simulations, which directly incorporate the molecular interactions underlying the 
wettability contrast on the flat substrate (see Fig. [21 for a snapshot). A number of different cases were considered (see 
Fig. which varied the amount of liquid present atop the stripe and the wettability of the exterior region. In order to 
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understand the numerical results, we have carried out a systematic linear stability analysis of the corresponding free 
boundary problem for the incompressible Navier-Stokes equation and find qualitative agreement. In addition to the 
full stability analysis we have presented a much simpler long-wavelength approximation for channel flow on straight 
chemical channels, which is able to capture some but not all of the qualitative features of the flow. 

The basic liquid configuration at rest, i.e., a non-volatile ridge of uniform cross-sectional shape, is known to be 
unstable due to a surface tension driven instability, similar to the Rayleigh-Platcau instability, if the contact angle is 
larger than 90° (see Fig.0. The MD simulat ions show that when the liquid is driven along the channel, the instability 
is enhanced by the driving if it is also present in the non-driven case (compare the times in Fig.[S]to those in Figs. Eland 
Fig. . However, statically stable liquid ridges are not destabilized by the flow. In the unstable case a periodic string 
of pearls appears growing in amplitude while propagating along the ridge. At later times, when nonlinear effects 
are manifest, the pearls develop distinct velocities and merge into pairs until only one of them survives. Despite 
substantial changes in amplitude and shape, the pearls remain confined to the chemical channel up to rather high 
accelerations. The instability has a beneficial side effect, in that the throughput is enhanced as compared to the 
homogeneously filled case (see Fig. EJ. For very high accelerations, however, the pearls are heavily deformed and in 
some cases lose contact with the substrate. A somewhat surprising result is that when the solid exterior to the wetting 
stripe is partially wetting rather than completely non- wetting, the qualitative behavior is unchanged (see Fig. |SJ) and 
the liquid remains atop the stripe. 

The wavelength of the pearling instability is reasonably well predicted by a full linear stability analysis of the 
Stokes free boundary problem (see Fig. ^| for the dispersion relations of the most unstable modes). For the acceler- 
ations studied in the MD simulations, the change in the wavelength of the most unstable mode as compared to the 
corresponding static case is too small to be detectable in MD simulations. However, the growth rate of the most 
unstable mode increases significantly upon forcing, in qualitative agreement with the MD simulations. The results are 
summarized in Tab. [I] A precise comparison between the MD results and the stability calculations is hampered by 
the thermal fluctuations present in the molecular simulation, which are exacerbated by the presence of a free surface 
which appears as a broad interfacial region (see Fig.^J. As a result, it is difficult to identify from the simulations the 
individual modes of the stability analysis, and their growth rates cannot be measured directly. 

The onset of the instability at contact angle 90° as well as the wavelength of the instability are reasonably well 
predicted by a simple long-wavelength expansion of the Stokes free boundary problem. Within this approximation, 
the only non-trivial part of the calculation is the numerical determination of the mobility (see Fig. Ill|) . However, 
since the most unstable wavelength is of the order of the ridge diameter, only qualitative results can be expected. In 
contrast to the MD results as well as to the full stability analysis, within this approximation the maximally unstable 
wavelength turns out to be independent of the acceleration (see Fig. ^] for the real part of the dispersion relation). 

The agreement between the MD simulation results and hydrodynamics (see Fig. [3| is less satisfactory than in other 
geometries without a free surface. The main reason seems to be the presence of thermal fluctuations of the liquid- 
vapor interface, i.e., capillary wave-like excitations, which broaden the liquid- vapor interface. The incorporation of 
thermal fluctuations has proven successful in relating MD simulations of a free liquid jet to the NS equations [22| . 
and more recently, have been included in a hydrodynamic thin film model p3l l24j . An extension of this technique to 
liquid ridges with large contact angles in chemical channels might improve the agreement with the simulations. 

Another possible source of disagreement between the MD simulations and the NS results is the neglect of the 
disjoining pressure in the hydrodynamic calculations. The issue arises because large portions of thin liquid films are 
exposed to the long-ranged substrate potentials and to the absence of liquid molecules outside the films [25ll2a| . and is 
expected to be important in the nanoscale flows studied here. In the simple case of a film of nearly constant thickness 
atop an infinite homogeneous substrate, the effect of the disjoining pressure is captured by adding a contribution 
—Atf/h 3 to the pressure, where Ah is the Hamaker constant, which is positive in a completely wetting case. If we 
repeat the long- wavelength analysis with this additional term added to the right hand side of Eq. I|29|) , the effect of 
the disjoining pressure is to modify the amplitude growth rate to 



7 ,ho 

WR = fl{— j .-1,2 

rj a a \ 



2(h 2 /a 2 - 1) 3A H ( a 



(h 2 /a 2 + l) 2 7 a 2 



+ arctan b t 
{kaf - (ka) 



Hence, the disjoining pressure has a stabilizing effect. It lowers the maximum growth rate and shifts the most unstable 

mode towards longer wavelength. In terms of the stability criterion, it leads to ho > a[l + (l + fr) 2 ] ^ 2 instead 

of ho > a. For realistic values of Ah 2(|, we get 3Ah/2"/ci 2 ~ 10 _1 . Such a difference will be masked by the 
uncertainties in MD simulations due to thermal noise and is undetectable at the present level of accuracy of our 
results. This calculation is only approximate, since the geometry in this paper is rather more complicated than a 
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nearly-uniform film on a homogeneous substrate. Other authors have extended this disjoining pressure analysis to 



simulations did not include a long-ranged disjoining pressure contribution per se, since the LJ interaction is cut off. 
The reason is that in MD simulations of flow at the 10 nm scale, the principal effect arises from the short-ranged 
interactions between liquid and solid atoms, and the long-ranged tail has a negligible effect (see, e.g., Ref. 



For wetting transitions however, this tail is quite important 25]. In particular on chemically structured surfaces, the 
influence of the laterally inhomogeneous substrate potential on the wetting film thickness has been discussed in detail 
p^. l30| and on chemical stripes morphological phase transitions have been predicted 0, . Recent experiments 
have confirmed the theoretical predictions on the structure of wetting films on a chemical stripe |33|. However, the 
thermodynamical ensemble plays a crucial role here. There are three types of ensembles which have to be distinguished 
due to their distinct characters. If the substrate is exposed to a large vapor reservoir one deals with a grand canonical 
ensemble as studied theoretically in Refs. [2jJ |3(1 I^H and experimentally in Ref. (3j|. For a volatile liquid enclosed 
into an isolating container one has a canonical ensemble; the MD simulations presented here correspond to this case. 
For a completely non- volatile liquid the liquid volume in the ridge is strictly conserved; this case has been studied in 
Refs. lalSl an d corresponds to the situation considered for the analytic hydrodynamic analyses in Sects. ITVl and Ivl 
The MD simulations are well described by the non-volatile Stokes dynamics, because the vapor pressure is very low 
and the vapor volume is small. Instabilities, such as the pearling instability, cannot be observed in a grand canonical 
setting (in accordance with the experimental observations in Ref. |33|), because they only occur in situations in 
which the pressure in the liquid ridge decreases with increasing volume. In this case, in a grand canonical setting, 
however, more and more molecules would be drawn out of the reservoir onto the ridge and it would grow without 
limit. Interestingly, the authors of Ref. [3^ find that the disjoining pressure does not influence the ridge shape above 
a ridge height of approximately 8 nm. This finding corroborates our conjecture given above, that the influence of the 
substrate potential on the instability studied in the present paper is small. 

In this work we have considered inertial driving, i.e., a force acting on every liquid atom in the system independent 
of its distance from the substrate surface. Although the physical mechanism for driving the liquid is qualitatively 
different, we expect electro-osmotically driven liquids to behave very similarly, in particular for low ion concentrations. 
If in such systems the Debye layer thickness is on the order of the ridge diameter, the ions are almost homogeneously 
distributed inside the ridge and the electrical force acts equally everywhere. This is in contrast to macroscopic systems, 
where the Debye layer is thin as compared to the system size and the electrical force only acts at the boundary of the 
system. 

Since the pearling instability increases the flow through a chemical channel (see Fig. [§J| , it might be an advantage 
to operate a microfluidic system in the unstable regime. Moreover, since the initially formed pearls have a relatively 
narrow size distribution, a chemical channel could be used as a nano-droplet dispenser. For longitudinal driving, 
the liquid stays on the chemical stripe even at high accelerations and when the instability is fully developed. For 
applications the stability with respect to non-longitudinal driving as well as the behavior at junctions and bends will 
be important. Full three-dimensional Stokes flow simulations, possibly taking into account the effect of long-ranged 
interactions and thermal fluctuations, should be used in order to interpret the results of MD simulations in such 
complicated geometries. 
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